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Abstract 

The article presents some aspects concerning the construction of a new thorn 
for the Cactus code, a complete 3-dimensional machinery for numerical relativity. 
This thorn is completely dedicated to numerical simulations in cosmology, that 
means it can provide evolutions of different cosmological models, mainly based on 
Friedman-Robertson- Walker metric. Some numerical results are presented, testing 
the convergence, stability and the applicability of the code. 

1 Introduction 

The core of modern cosmology is the theory of General Relativity, a 4-dimensional 

theory involving one dimension of time and three of space, having as field equations the 
Einstein equations : 

1 8?rG 
Rij - 2 gi i R + - ~^4~ Ti i CO 

where A is the cosmological constant, R^ the Ricci tensor, R the Ricci scalar, the 
space-time metric, Ty the stress-energy tensor, G the gravitational constant, c the speed 
of light and i,j = 0,1,2,3. Numerical Relativity (NR) is concerned with the study 
of numerical solutions of the Einstein's equations for the gravitational field (||), and 
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when we apply these equations to cosmological models we have numerical cosmology as 
a branch of NR. Einstein equations are an extremely complicated system of coupled, 
non-linear, partial differential equations and solving them numerically makes enormous 
demands on the processing power and memory of a computer. The Cactus code (||) 
has been mainly designed as a computational toolkit (freely available for the scientific 
community) for simulating different systems of partial differential equations, as Einstein 
equations are The code can be used to simulate fully 3-dimensional systems with strong 
gravitational fields: collapsing gravitational waves, colliding black-holes, neutron stars, 
and other violent astrophysical processes generating gravitational waves. 

As a branch of NR, numerical cosmology can use the Cactus code for his purposes. Al- 
though Einstein field equations reduce to a small number of equations in most cosmological 
models (0]-|§), trying to simulate more complex models (as models with cosmological 
constant or having matter fields coupled with gravity) it becomes difficult to compose 
special numerical codes for every specific case. Thus our purpose was to develop a piece 
of code to be used inside a general numerical machinery for solving Einstein equations, 
as the Cactus code is. This article is a report on how we extended the Cactus code 
for numerical cosmology. In this purpose we developed a new thorn of the Cactus code 
( "thorns" are called the computer code packages embedded in Cactus code for some spe- 
cific application). The starting point of this new thorn, called "Cosmo" was an existent 
thorn of the Cactus code (the "Exact" thorn) previously reported in (||) for testing the 
Cactus code on exact solutions of the Einstein equations. 

We investigated some cosmological models being exact solutions of the Einstein equa- 
tions in our previous article (f|). Here the scale factor of the universe (R(t)) and thus the 
components of the metric tensor are known analytic function of time, given in the code 
at the initial time as initial data. But in our new Cosmo thorn, the main point is that we 
have to use the general Robert son- Walker metric, having the scale factor of the universe 
(R(t)) as a unknown function to be evolved during the time evolution of the Einstein 
equations through the Cactus code. Thus the only job to be done by the Cosmo thorn is 
to provide good initial data (that means at the initial time to considered the actual time 
of the universe, for example) and then let the Cactus code to solve the Einstein equations 
through the time. In addition it is necessary to prescribe convenient boundary values 
for the functions and variables involved. This is necessary because the Cactus code has 
a limited boundary methods implemented, not appropriate for our numerical cosmology 
purposes. 

The article is organized as follows : the next section 2, presents some general facts 
about the way we built the Cosmo thorn, how we solved the problems exposed above 
and the necessary theoretical cosmology background. Section 3 is dedicated to a matter 
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dominated universe model based on Robert son- Walker metric (10]), with pressure- less 
matter - "dust" and without cosmological constant. Models with cosmo logical constant 
are presented in the next section 4 where we concentrated again only on models having 
matter as a pressure- less dust. A wide larger class of possible initial data is possible 
including the so called "coasting cosmologies", models not having as a starting point a 
Big Bang. In all cases studied and reported in this article we investigated the convergence 
and the stability of the Cactus code, in the same way we presented in detail in ||. 

As a major conclusion we have to point out the good second order convergent and 
stable behavior of the code, even in the cases we studied both forward and backward 
time evolutions. Future perspective are open for developing this thorn for inflationary 
models. It is in our view to study those models having one or more scalar field coupled 
with gravity to control the behavior of the model and to simulate an accelerated late time 
evolution, actually called "cosmic acceleration" being well proved experimentally in the 
last years. Thus together with other similar thorns being in preparation at AEI, we shall 
have a new class of applications for the Cactus code, a major instrument for new studies 
in numerical cosmology and in cosmology, in general. 

Being included in the list of thorns of the Cactus code, actually the files for our thorn 
are available by request, directly at the author's e-mail address and in the near future it 
will be included in the CVS repository of the Cactus code (see |J) for free download. 

The simulations were performed mainly on single processor machines, using both the 
AEI computer network (SGI or Dec machines with UNIX operating system) and a Pen- 
tium III machine with a UNIX FreeBSD operating system at the West University of 
Timi§oara. Some of the simulations, with similar results, were also done on the Origin 
2000 supercomputer at the AEI. 

Through this article and in the Cactus code we shall use geometrical units with G = 
c= 1. 



2 The Cactus code, the Robertson- Walker metric 
and the Cosmo Thorn 

As we mentioned earlier, in modern cosmology we are using the Robertson- Walker metric 
(RW), namely 



d s 2 = - c 2 dt 2 + R(t) 



dr 2 
1 — kr 2 



+ r 2 (d6 2 + sin 2 6 dcj) 2 ) 



(2) 
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or in isotropic coordinates 



ds 2 = -c 2 df 

where 



R(ty 



l + fr 



r' 



dr' 2 + r'\d0 2 + sin 2 8 dtf) (3) 



2r 



1 + VI - kr 2 



as a generic metric for describing the dynamics of the universe through the Einstein 
equations ([I]). Here A; is a constant with arbitrary value, positive (for closed universes), 
negative (for open universe) and zero for flat universes. Usually, this constant is taken 1, 

— 1 or respectively. When we consider all the matter in the universe as a perfect fluid, 
the stress-energy tensor may be given as 

Tij = (e + p)uiUj + pgij (4) 

and is necessary to prescribe a state equation, i.e. a relation between the pressure p and 
the density of the universe p. For matter dominated universes usually the pressure is 
p = and this will be the case we concentrated in our simulations. 

As our main purpose was to build a new thorn for the Cactus code, generically named 
from now one "The Cosmo thorn" for evolving the RW metric above in all three cases for 
k = —1,0 and +1. The main problem was to to "convince" Cactus to evolve, without 
having the time behavior of R(t) function (actually, R(t) is a function of the thorn). 
Although we used the same strategy as in the above mentioned Exact thorn (see f6j) 
here, not having an analytic function of time for R(t) the solution was to provide Cactus 
code all initial data he needs, namely R(t) and R(t) = dR(t)/dt at the initial time to we 
choose ! As a result, new routines were added for calculating the extrinsic curvature J51 
components at the initial time, namely : 

K aP = -2R{t) 2 ^-5 aP = -2R(t) 2 H(t)6 aP (5) 
R{t) 

where H(t) is the Hubble constant function (!). Note that we have here as lapse 
function, N — 1 and a shift vector as N a — 0, a, (3 — 1, 2, 3 as always in RW cosmology. 

Another item we concentrated on was the boundary problem : - Cactus code and the 
thorns which solve numerically the Einstein equations (namely ADM_BSSN and ADM 

- see p|) have not implemented proper boundary conditions for the RW metric. In our 
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previous article on Exact thorn we used "flat" boundary condition but this is not 
appropiate here ! To solve this problem we used a simple hint : R(t) and H(t) are 
functions of time only, thus they are constant on the numerical grid at one time, having 
the same value inside the grid and on the boundary. As a result, we added a new boundary 
method through a new value of the parameter bound of the ADM_BSSN thorn (see || 
and ||), called external which activates some new routines we composed in this purpose, 
having as a generic name recover. In these routines the spatial metric components g a @ 
and the extrinsic curvature components K a p on the boundaries are "recovered" from the 
values of the R(t) and H(t) on the interior points by calculating them using equations 
(0) and (||). Also these new boundary conditions are necessary during the intermediate 
steps of integration used by ADM_BSSN thorn in the so-called "IterativeCN.F" routine 
(where a three steps Cranck-Nicholson integration method is applied). For this purpose 
we slightly modified the "IterativeCN.F" routine for "injecting" our specific "recover" 
boundary values. Actually this was possible through a facility of the Cactus code which 
makes possible to call different routines from different thorns using special "include" 
commands (see ||). 



3 Friedman-Robertson- Walker (FRW) cosmologies 

By "FRW cosmologies" we mean those cosmological models based on the above RW metric 
and being solutions of the Einstein equations (0) without cosmological constant (A = 0) - 
see for example (J7|) and (f8|). Restricting ourselves only to matter dominated universes 
(i.e. p(t) = 0) as we mentioned earlier, and solving the conservation law of the perfect 
fluid which emulates the matter, we have for the stress-energy tensor components 

c 2 C 

T a/3 = and T 00 = — ^ (6) 

where C is a constant depending on the the mass density of the universe at the initial 
time to- Solving the Einstein equations for this constant and at the initial time to where 
R(t ) = Ro and H(to) = H , we have 

C = n oPc Rl = 2qo Pc Rl = goc 2 4 ^; )3 (7) 
where go is the deceleration parameter at the initial time, i.e. 



Q(t) 



H{t) 2 R(t) 
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and Qq is the so-called density factor, defined as 



(9) 



Po being the initial density of the universe and p c the critical density (i.e. the density of 
a flat three-dimensional universe where k = 0). It is proven (||) that we have : 



flo = 2q and p c 



3H[ 
8nG 



(10) 



In order to provide initial data for the Cactus code we need to solve the Einstein equations 
at the initial time for the initial Hubble constant Hq (to calculate the values of the extrinsic 
curvature through equation (||)) and for the initial scale factor R . It is known that we 
have here three solutions, depending on the geometric factor k = 0,1, —1. Thus, for k = 
(flat three-dimensional spacetime) where 



we have (see ?? eq. 4.37) 



and Qq = 1 



H = — 
3t 



and we shall take Rq = 1, as a value for calculating, for now one the relative scale factor, 
i.e. R(t)/RQ which will be our numerical output. For the case of closed universe, k = 1 
we have 



go > 



and Qq > 1 



In this case (as in the next one for k = — 1) we do not have an analytic solution for the 
scale factor, but for the initial time to we can easily obtain that (?? eq. 4.51) : 



Hn 



qo 



(2g - If' 2 
and the initial scale factor is : 



cos 



qo 



qo 



V2q^T 

qo 



Ra 



Finally, for the open universes (k 



H oy /2q ~~T 
where we have 



1 

T 



'12) 



(13) 



qo < 



and Qq < 1 
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As it is shown (|| eq. 4.64), here we have : 

" 1 



a * 



(1 - 2q„)W 
and the initial scale factor is : 



VI - 2g _ ln / 1 - g + yi - 2g 
9o V 9o 



(14) 



* = iWT^S (15) 

With these results reminded, we can now proceed to explain how we produced our sim- 
ulations. Taking as initial data the initial time (to) in billion years, being an internal 
parameter of the Cactus code, namely CCTK_intitial_time, we then have all the vari- 
ables (in our geometrical units with c = G = 1, of course) rescaled in billion years. 
We also have the initial value of the deceleration parameter, go denoted in the code as 
robson_q. As output we have the scale factor, more precisely the relative scale factor 
R(t)/R (denoted in the code as raza ) and the Hubble parameter function H(t) (denoted 
in the code as hubble. For evolving these two functions we produced two internal rou- 
tines of the thorn, which mainly calculates their values from the general output values 
provided by the Cactus code and the ADM_BSSN thorn (as the extrinsic curvature and 
the three-dimensional metric components). 

Here we shall present the results we obtained on evolving a matter dominated universe 
(i.e pressure-less matter, p — 0) based on the RW metric. Some explanations are necessary 
in order to clarity our figures ( see also ([§[))■ Here and in the next figures, "normalized" 
means the L2 norm of the respective function calculated using all the values of the function 
on the computational grid at one time and the respective output Cactus files are denoted 
with _nm2 extension. 

Cactus code is evolving the Einstein equations using the 3+1 decomposition of space- 
time (both in ADM or BSSN evolution methods -see (||) and (|[L0|)). Thus the Einstein 
equations can be split in two groups : the dynamic equations (for the time derivatives of 
the three-dimensional metric and extrinsic curvature) and the constraint equations (the 
Hamiltonian constraint and the Momentum constraint) - see |IJ and 0. The constraint 
equations are satisfied during all time evolution of the system. Thus one of the main tests 
on the convergence in the Cactus code is given by the time behavior of the Hamiltonian 
constraint (an output of the Cactus code through a thorn called ADMConstraints, namely 
ham). In some of the next figures (here and in the next sections) we show the convergence 
of the L2 norm of the Hamiltonian constraint for different number of iterations, using two 
different resolutions on the grid, one with 14 3 and the second with 28 3 points (both grids 
cover the same region of spacetime, so the grid with more points has a smaller value of the 
finite difference interval Ax). Notice that the Hamiltonian constraint for a true solution 
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Figure 1: Convergence test for FRW flat matter dominated universe (k=0, p=0) for 
100/200 iterations (left side) and for FRW closed matter dominated universe (k=l, p=0) 
1000/2000 iterations (right side) 
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of the Einstein equations should be equal to zero. Finite differencing errors imply that 
the numerical solution will have a non- vanishing value of the Hamiltonian constraint. For 
a consistent finite difference approximation of the Einstein equations, we should expect 
the Hamiltonian constraint to approach zero as the resolution is increased. For a second 
order approximation, the value of the Hamiltonian constraint should go down by a factor 
of four when the resolution is doubled. This was the method to test the convergence of 
the code in all examples presented in this article. We have obtained good second order 
convergence in all examples as is shown in the next figures. 

First of our figures (|l] and ||]) are pointing out the convergence of the code for different 
examples of FRW cosmologies. The time is scaled in billion years (as for all figures in 
this section) starting with an initial time of 20 billion years. As it is obvious from these 
figures a good second order convergence is visible, even for long time evolutions. For the 
case of the closed universe (k = 1) we had to restrict our simulations only till around 
9000 iterations (for a grid with 14 3 points) because in this case the universe is evolving 
to a Big Crunch (opposite to the Big Bang) as it is pointed out in the right panel of the 
next figure @. Here the scale factor is evolving to a maximum having the value 

_ 2g _c_ 
max (2g - 1)3/2 ■ Hq 

twice of the actual value Ro (we used here go = 1)- For spatially flat universes we used, 
of course, go = 0.5 and for open universes we had go = 0.25. 

Figure (|3|) contains also in his left panel the evolution of the Hubble parameter function 
in time. We plotted here, for convenience the modulus of this function being calculated 
through the L2 norm as is output by the code. 

The same comment is valuable for the left panel of the figure @ where we can point 
out also the coincidence of the minimum zero value of the Hubble parameter with the 
maximum of the scale factor from the previous figure @. The right panel of (§) is 
dedicated to the simultaneous plot of the three cases, pointing out the differences between 
the geometries of the three models. 

Last of our figures from this section, (^) contains the plots of the relative scale factor 
time evolution for different closed matter dominated FRW universes, having different 
deceleration parameter (i.e. different density parameters). It is visible here that as the 
deceleration parameter (density parameter) is increasing so is decreasing the time-life of 
the closed universe between a Big-Bang and a Big-Crunch. 
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Figure 2: Convergence test for FRW open matter dominate universe (k=-l, p=0) for 
1000/2000 iterations (left side) and for 9000 iterations (right panel) 

4 FRW cosmologies with cosmological constant 

Here we shall deal with cosmologies based on the above RW metric but having a non- 
vanishing cosmological constant. In spite of the failure of the the initial tentative of 
Einstein to use the cosmological constant as a trick to fit general relativity with, at 
that time presumed, static universe, the cosmological constant is again used in modern 
cosmology in order to establish more realistic models of the universe, as the inflationary 
models and, recently, the so called cosmic acceleration. Early attempts (see (@)) on 
using numerical routines for solving the Einstein equations for these models were reported. 
Because our goal is to implement in the Cactus code a thorn for doing similar simulations, 
but on a larger scale and fully three-dimensional we shall refer to these previous results 
only for comparing our results. Meanwhile we shall use here the same notations and 
scaling of the parameters as in (||) for prescribing the initial data sets for the Cactus 
code. 

This time I will proceed in a slightly different way than in the previous section, scaling 
all our variables and parameters in terms of the initial value of the Hubble parameter 
function H , which is not completely free in our universe, having been measured to within 
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Open Friedmann universe (k=-1 , p=0) Closed Friedmann universe (k=1, p=0) 

Time evolution of the Hubble function Time evolution of the scale factor (9000 iterations) 




20 25 30 35 40 50 100 150 200 

time (billion years) time (billion years) 



Figure 3: Time evolution of the Hubble function for FRW open matter dominated universe 
(k=-l, p=0) 1000/2000 iterations (left panel) and of the relative scale factor for for FRW 
closed matter dominated universe (k=+l, p=0) 9000 iterations (right panel) 
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Closed Friedmann universe (k=1 , p=0) 

Time evolution of the Hubble function (L2 norm, 9000 iterations) 

0.03 



Friedmann matter dominated universes 

Time evolution of the relative scale factor 




100 

time (billion years) 



1 00 150 200 

time (billion years) 



Figure 4: Time evolution of the Hubble function (L2 norm) for FRW closed matter 
dominated universe (k=+l, p=0) 9000 iterations (left panel) and the time evolution of 
the relative scale factor for FRW matter dominated universes (p=0, k=+l,0,-l) 9000 
iterations (right panel) 
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Figure 5: Time evolution of the relative scale factor for for different closed matter domi- 
nated FRW universes (k=+l, p=0) 9000 iterations 

a factor of 2. Thus for the cosmological constant A for instance, we use X/Hq as a 
basic parameter, denoted in the code with robsonJ. Of course, for dimensional scaling 
reasons we prescribed at the initial time Hq = 1. Thus, for example, at the output we have 
the relative value of the Hubble function H(t)/H~Q. Meanwhile we need a parameter to 
describe the density of the universe, and this time we used fl instead of the deceleration 
parameter (remark that here we have no more the relation Q = 2q ). We have also the 
Gaussian curvature K(t) defined as 

which will be, of course identical to the geometric factor k at the initial time when we 
take the initial scale factor Rq = 1 as we done in our simulations. Following the line 
indicated in we can define the density factor as 

87rG pit) 8ttG p 

il[t) = — — and at the initial time iZ = — — 775 (16) 

3 H(t) 3 Hq 

where po is the initial density of the matter in the universe. Now solving the Einstein 
equations at the initial time and after some algebraic manipulations we obtain (for c = 
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G = l,Ro = l): 

K = k = Hl(n-l+ 1 --^ (17) 

and for the stress-energy tensor (where we included also the term containing the cosmo- 
logical factor, as we did in (0)) : 

TiJ = ^ " A ^ Where T °°° = ^G n ° H2 °i§f ; T ^ = ° V "' ^ = 2 ' 3 (18) 

Now our numerical results are pointed out in the next figures. First we show, in figure (§) 
the time evolution of the scale factor for some cosmological models, having the density 
factor Q being 0.1 (left panel) and 1.0 (right panel). The models are differentiated by 
their respective cosmological constant (more precisely the factor X/Hq). For the time 
scale we used here Hubble time units, i.e. the time is expressed as tjtu where tjj is the 
Hubble time : tu = Hq . In these simulations, as in the next ones we started from the 
actual time (taken here to — 0) and we investigated the time behavior of the model in 
both directions of time, forward, as well as backward in time. 

Next figure (|7|) shows, in his left panel some special cases of FRW models with cosmo- 
logical constant, where the scale factor shows off a different evolution : namely it starts 
somewhere back in time with a certain finite value, it decreases in time, then it has a 
period (longer or shorter) of constant value, evolving finally as an expanding model (of 
De-Sitter type, for example). These models, having a period of constant value of the scale 
factor, are sometimes called "coasting cosmologies" and are investigated as an alternative 
(realistic or not) to the standard Big-Bang model, not being generated from a Big-Bang, 
as can be seen from our simulations. The right panel of the figure (|7|) contains the time 
evolution of the Hubble constant function (the L2 norm, as it is output by the Cactus 
code, otherwise the values after the null minimum value being negative ones) for a model 
with Q = 0.1 and \/H$ = 5. 

We must point out that we performed here also several convergence tests, as in the 
previous section, obtaining again good second order convergence. As an example we shall 
present in the last figure @ the convergence test for 14 3 and 28 3 points on the grid for a 
model having again Q — 0.1 and X/Hq = 5 for an evolution forward in time (left panel) 
and backward in time (right panel). 
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Friedmann models with cosmological constant Friedmann models with cosmological constant 



Density parameter O.I Density factor 1 




Cosmic time (in Hubble time units) Cosmic time (in Hubble time units) 



Figure 6: Time evolution of the scale factor for FRW models with cosmological con- 
stant. The left panel represent four different models for a density factor Q = 0.1 and 
A/Hq = {—3,0,3,6} and the right panel represent models for a density factor Q = 1.0 
and \/Hl = {-1,-0.5,9} 
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Friedmann models with cosmological constant Friedmann models with cosmological constant 



Coasting cosmologies with density factor 0.1 Density parameter 0.1 




Cosmic time (in Hubble time units) Cosmic time (in Hubble time units) 

Figure 7: "Coasting" FRW cosmologies for Q = 0.1 and X/Hq = 4.94,5,6 (left panel) 
and the time evolution of the Hubble constant function for a model with fl = 0.1 and 
X/Hq = 5 (right panel) 
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Figure 8: Convergence test for 14 3 and 28 3 points on the grid for a model having again 
Q = 0.1 and \/Hq = 5 for an evolution forward in time (left panel) and backward in time 
(right panel) and for 100/200 iterations 
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